# Replication materials for 
# Peisker, Hoffmann, Muttarak (2025)
# Climate news mediates extreme weather effects on climate change concern

packages <- c(
  "tidyverse", "broom", "data.table", "parallel",
  "mediation", "texreg", "xtable", 
  "scico", "patchwork", "ggplotify"
)
lapply(packages, library, character.only = TRUE)
lapply(packages, citation)

source("R Files/functions.R")
theme_set(theme_bw())
setDTthreads(parallel::detectCores())
options("mc.cores" = detectCores())

load("Data/eb_climate_reg.RData")

#### baseline specification ####
#### mediator models ####
yv <- "emm_climate_w_nat_4w"
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
)
fml_list_med <- lapply(xv, function(x){
    f_lm(yvar = yv, xvars = x)
  })
tab_mediator <- list(
  lm(fml_list_med[[1]], data = eb_climate_outcome_c),
  lm(fml_list_med[[2]], data = eb_climate_outcome_c),
  lm(fml_list_med[[3]], data = eb_climate_outcome_g)
)
screenreg(tab_mediator)

#### outcome models ####
yv <- "cc_eu_concern_w"
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
)
fml_list_out <- lapply(xv, function(x){
  f_lm(yvar = yv, xvars = x)
})
tab_outcome <- list(
  lm(fml_list_out[[1]], data = eb_climate_outcome_c),
  lm(fml_list_out[[2]], data = eb_climate_outcome_c),
  lm(fml_list_out[[3]], data = eb_climate_outcome_c),
  lm(fml_list_out[[4]], data = eb_climate_outcome_g)
)
tab_export <- c(tab_mediator, tab_outcome)
my_screenreg(tab_export)

#### mediation analysis ####
tab_mediate_pars <- mediate(
  tab_mediator[[1]], tab_outcome[[2]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_pars)

tab_mediate_control <- mediate(
  tab_mediator[[2]], tab_outcome[[3]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_control)

tab_mediate_control2 <- mediate(
  tab_mediator[[3]], tab_outcome[[4]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_control2)

# save(tab_mediate_pars, tab_mediate_control, tab_mediate_control2, file = "Data/tab_mediate_rn.RData")
# load("Data/tab_mediate_rn.RData")
pars_tidy <- 
  tab_mediate_pars %>%
  tidy_med() %>% 
  mutate(
    Specification = "parsimonious"
  )  
control_tidy <- 
  tab_mediate_control %>%
  tidy_med() %>% 
  mutate(
    Specification = "with controls"
  )  
control2_tidy <- 
  tab_mediate_control2 %>%
  tidy_med() %>% 
  mutate(
    Specification = "with controls 2"
  )  

#### sensitivity ####
tab_mediate_pars_sens <- 
  medsens(tab_mediate_pars, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_pars_sens)
tab_mediate_control_sens <- 
  medsens(tab_mediate_control, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_control_sens)
tab_mediate_control2_sens <- 
  medsens(tab_mediate_control2, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_control2_sens)
# save(tab_mediate_pars_sens, tab_mediate_control_sens, file = "Data/tab_mediate_sens.RData")

pdf(file = "Figures/medsens_rn.pdf", width = 10, height = 5)
par(mfrow=c(1,2)) 
plot(tab_mediate_pars_sens, sens.par = "rho", xlim = c(-0.4,0.4), ylim = c(-0.2,0.2))
plot(tab_mediate_pars_sens, sens.par = "R2", sign.prod = "positive", r.type = "residual", xlim = c(0,0.5), ylim = c(0,0.5))#
dev.off()

#### make table ####
tab_med_tidy <- 
  pars_tidy %>% 
  bind_rows(control_tidy) %>% 
  bind_rows(control2_tidy) %>% 
  mutate(
    estimate = paste0(round(estimate, 3), "^{", stars, "}"),
    ci_95 = paste0("(", round(ci_low, 2), ", ", round(ci_high, 2), ")")
  )  %>% 
  pivot_longer(cols = c(estimate, ci_95), names_to = "name")

tab_med_tidy2 <- 
  with(tab_med_tidy,
       tibble(
         term = c("Mediator model", term[1:8]),
         col1 = "", col2 = "", col3 = "", col4 = "",
         col5 = c("(1)", value[Specification == "parsimonious"]),
         col6 = c("(2)", value[Specification == "with controls"]),
         col7 = c("(3)", value[Specification == "with controls 2"])
       ))
tab_med_tidy2$term[seq(3, 9, 2)] <- ""
tab_med_tidy2

tab_sens <- cbind(
  term = c("$\\rho$", "$R^{*2}_{M}R^{*2}_{Y}$", "$\\tilde{R}^2_{M}\\tilde{R}^2_{Y}$"),
  col1 = "", col2 = "", col3 = "", col4 = "",
  col5 = tidy_medsens(tab_mediate_pars_sens),
  col6 = tidy_medsens(tab_mediate_control_sens),
  col7 = tidy_medsens(tab_mediate_control2_sens)
  ) 
tab_sens 
tab_med_final <- 
  tab_med_tidy2 %>% 
  rbind(c("Sensitivity analysis", rep("", 7))) %>% 
  rbind(tab_sens)
tab_med_final
print(
  xtable(tab_med_final), #
  file = "Tables/tab_med_tidy.tex",
  include.rownames = FALSE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function	= identity
)

#### with spatio-temporal polynomials ####
#### mediator models ####
cpol_vars <- grep("_lon", names(eb_climate_outcome_c_sp), value = T) 
yv <- "emm_climate_w_nat_4w"
xv <- paste0(c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
), " + ", paste0(cpol_vars, collapse = " + "))

fml_list_med <- lapply(xv, function(x){
  f_lm(yvar = yv, xvars = x)
})
tab_mediator_sp <- list(
  lm(fml_list_med[[1]], data = eb_climate_outcome_c_sp),
  lm(fml_list_med[[2]], data = eb_climate_outcome_c_sp),
  lm(fml_list_med[[3]], data = eb_climate_outcome_g_sp)
)
screenreg(tab_mediator_sp)

#### outcome models ####
yv <- "cc_eu_concern_w"
xv <- paste0(c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
), " + ", paste0(cpol_vars, collapse = " + "))
fml_list_out <- lapply(xv, function(x){
  f_lm(yvar = yv, xvars = x)
})
tab_outcome_sp <- list(
  lm(fml_list_out[[1]], data = eb_climate_outcome_c_sp),
  lm(fml_list_out[[2]], data = eb_climate_outcome_c_sp),
  lm(fml_list_out[[3]], data = eb_climate_outcome_c_sp),
  lm(fml_list_out[[4]], data = eb_climate_outcome_g_sp)
)
tab_export <- c(tab_mediator_sp, tab_outcome_sp)
my_screenreg(tab_export, omit.coef = "_lon")

#### mediation analysis ####
tab_mediate_sp_pars <- mediate(
  tab_mediator_sp[[1]], tab_outcome_sp[[2]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_sp_pars)

tab_mediate_sp_control <- mediate(
  tab_mediator_sp[[2]], tab_outcome_sp[[3]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_sp_control)

tab_mediate_sp_control2 <- mediate(
  tab_mediator_sp[[3]], tab_outcome_sp[[4]], 
  treat = "tmean_nf_pos1_4w",
  mediator = "emm_climate_w_nat_4w",
  boot = TRUE, boot.ci.type = "perc",
  sims = 1e4
)
summary(tab_mediate_sp_control2)

# save(tab_mediate_sp_pars, tab_mediate_sp_control, tab_mediate_sp_control2, file = "Data/tab_mediate_sp_rn.RData")
# load("Data/tab_mediate_sp_rn.RData")
sp_pars_tidy <- 
  tab_mediate_sp_pars %>%
  tidy_med() %>% 
  mutate(
    Specification = "parsimonious"
  )  
sp_control_tidy <- 
  tab_mediate_sp_control %>%
  tidy_med() %>% 
  mutate(
    Specification = "with controls"
  )  
sp_control2_tidy <- 
  tab_mediate_sp_control2 %>%
  tidy_med() %>% 
  mutate(
    Specification = "with controls 2"
  )  

#### sensitivity ####
tab_mediate_sp_pars_sens <- 
  medsens(tab_mediate_sp_pars, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_sp_pars_sens)
tab_mediate_control_sp_sens <- 
  medsens(tab_mediate_sp_control, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_control_sp_sens)
tab_mediate_control2_sp_sens <- 
  medsens(tab_mediate_sp_control2, rho.by = 0.1, effect.type = "indirect", sims = 1e4, eps = 0.001)
summary(tab_mediate_control2_sp_sens)
# save(tab_mediate_pars_sens, tab_mediate_control_sens, file = "Data/tab_mediate_sens.RData")

pdf(file = "Figures/medsens_sp_rn.pdf", width = 10, height = 5)
par(mfrow=c(1,2)) 
plot(tab_mediate_sp_pars_sens, sens.par = "rho", xlim = c(-0.4,0.4), ylim = c(-0.2,0.2))
plot(tab_mediate_sp_pars_sens, sens.par = "R2", sign.prod = "positive", r.type = "residual", xlim = c(0,0.5), ylim = c(0,0.5))#
dev.off()

#### make table ####
tab_med_sp_tidy <- 
  sp_pars_tidy %>% 
  bind_rows(sp_control_tidy) %>% 
  bind_rows(sp_control2_tidy) %>% 
  mutate(
    estimate = paste0(round(estimate, 3), "^{", stars, "}"),
    ci_95 = paste0("(", round(ci_low, 2), ", ", round(ci_high, 2), ")")
  )  %>% 
  pivot_longer(cols = c(estimate, ci_95), names_to = "name")

tab_med_sp_tidy2 <- 
  with(tab_med_sp_tidy,
       tibble(
         term = c("Mediator model", term[1:8]),
         col1 = "", col2 = "", col3 = "", col4 = "",
         col5 = c("(1)", value[Specification == "parsimonious"]),
         col6 = c("(2)", value[Specification == "with controls"]),
         col7 = c("(3)", value[Specification == "with controls 2"])
       ))
tab_med_sp_tidy2$term[seq(3, 9, 2)] <- ""
tab_med_sp_tidy2

tab_sens_sp <- cbind(
  term = c("$\\rho$", "$R^{*2}_{M}R^{*2}_{Y}$", "$\\tilde{R}^2_{M}\\tilde{R}^2_{Y}$"),
  col1 = "", col2 = "", col3 = "", col4 = "",
  col5 = tidy_medsens(tab_mediate_sp_pars_sens),
  col6 = tidy_medsens(tab_mediate_control_sp_sens),
  col7 = tidy_medsens(tab_mediate_control2_sp_sens)
) 
tab_sens_sp
tab_med_sp_final <- 
  tab_med_sp_tidy2 %>% 
  rbind(c("Sensitivity analysis", rep("", 7))) %>% 
  rbind(tab_sens_sp)
tab_med_sp_final
print(
  xtable(tab_med_sp_final), #
  file = "Tables/tab_med_tidy_sp.tex",
  include.rownames = FALSE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function	= identity
)

